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Abstract 

Power  spectral  estimation  provides  one  approach  to 
direction-finding.  This  approach  readily  generalizes 
to  produce  a  collection  of  newer  direction-finding  algo¬ 
rithms.  Estimation  of  the  bispectrum  yields  a  bispec- 
tral  direction  finder.  Estimates  of  time-frequency  dis¬ 
tributions  produce  Wigner-Ville  and  Gabor  direction¬ 
finders.  Some  types  of  non-stationary  time  series  ad¬ 
mit  spectral  estimators  which  can  be  used  to  localize 
a  source.  Chaotic  signals  can  also  be  localized  using 
recently  developed  parameter  estimators. 

1  Introduction 

In  1982,  Don  Johnson  wrote  a  much-cited  paper  en¬ 
titled  The  Application  of  Spectral  Estimation  Methods 
to  Bearing  Estimation  Problems  (10).  In  the  interven¬ 
ing  years,  considerable  development  has  taken  place 
in  time-frequency  distributions,  higher-order  statis¬ 
tics,  and  estimators  for  non-stationary  random  pro¬ 
cess.  This  paper  extends  the  theme  developed  by 
Johnson  to  encompass  some  of  these  newer  results. 
Following  a  review  of  conventional  direction-finding 
(DF),  section  3  shows  how  estimators  of  the  bispec¬ 
trum  can  yield  DF  algorithms.  DF  using  the  Gabor 
Transform  and  Wigner-Ville  Distribution  (WVD)  are 
covered  in  sections  4  and  5.  Spectral  estimators  for 
time  series  which  are  periodically  correlated  (cyclo¬ 
stationary)  ako  yield  a  DF  algorithm  (section  61  as 
do  the  estimators  for  a  fractrd-like  time  series  (sec¬ 
tion  7). 

Vectors  are  denoted  boldface  x.  The  conjugate 
transpose  is  denoted  as  .  The  normal  or  Gaussian 
distribution  with  mean  fi  and  variance  is  denoted 
//{n,cr^).  If  {x(<)}  is  a  zero-mean,  vector-valued  ran¬ 
dom  process,  the  covariance  function  is 


/2x(ti,<2)  =  £^[x(<i)x(t2)"], 


and  the  cross  spectral  density  (CSD)  matrix  is: 


Px(h,h) 


Table  1:  Principal  Notation 


t 

time 

f 

frequency 

c 

speed  of  sound 

E 

expectation  operator 

u 

unit  (look)  vector 

Xm(i) 

mth  sensor  output 

H 

II 

array  “snapshots” 

r  m 

vector  to  mth  sensor 

a(/,u)  =  [exp(i2T/r^u/c)] 

steering  vector 

2  Conventional  DF 

A  quick  review  of  conventional  DF  is  undertaken 
in  this  section  to  provide  the  nomenclature  for  the 
subsequent  DF  algorithms.  Assume  a  single  source 
is  transmitting  a  signal  s(<)  in  a  lossless,  plane-wave 
environment  in  the  direction  — uq.  The  signal  s{t)  is 
received  by  M  omni-directional  sensors  fixed  at  posi¬ 
tions  r,„.  Each  sensor  receives  a  time-delayed  copy  of 
s(t).  The  time  delay  function  associated  with  the  mth 
sensor  is  r„,(uo)  =  r^Uo/c.  If  additive  noise  {z„{t)) 
is  present  at  the  mth  sensor,  then  the  noise-corrupted 
signal  at  the  m-th  sensor  is  determined  as  (10): 

*m(0  =  «(<  +  T-m(Uo))  + Zm(<)-  (D 

The  classical  DF  output  x(t,  u)  using  look  vector  u  is 
determined  by  weights  {wm)  and  time  delays  (rm(u)) 
as 

M 

x(t,u)  =  ^  VJmXmit  “  ^"^(u)).  (2) 

m  =  1 

The  fundamental  point  of  Johnson’s  paper  [10]  is  that 
the  estimates  of  the  power  spectrum  of  {x(t,ujl  yield 
bearing  estimates.  To  see  this,  make  the  following 
standard  assumptions  regarding  the  stationarity  of  the 
signal  and  the  noise  process  z(l)  [rm(0]' 


Assumption  2.1  Single  source,  uncorrelated  noise: 

1.  the  plane-wave  signal  s{t)  is  a  zero-mean,  wide-sense 
stationary  (WSS)  random  process. 

2.  the  noise  z{t)  is  a  zero-mean,  wide-sense  stationary 
random  process. 

S.  the  signal  j(t)  and  the  sensor  noise  are  uncorrelated. 

If  s{t)  WSS,  it  admits  a  mean-square  Fourier  rep¬ 
resentation 

s{t)  =  r 

J  — oo 

where  dS{f)  is  an  orthogonally  scattered  measure  [12]. 
The  power  spectrum  P,  {f)  of  the  signal  is  determined 
by  dP,{f)  =  F?[|dS(/)p].  Likewise,  the  WSS  noise 
vector  z(t)  admits  Fourier  representation  with  an  or¬ 
thogonally  scattered  vector  measure  dZ(/)  and  its 
CSD  matrix  is  given  by  dPz{f)  =  E[dZ(f)d7i{f)^]. 
Under  the  assumptions  2.1,  it  may  be  shown  that  the 
array  output  {x(t)J  admits  an  orthogonally  scattered 
vector  measure  [lOJ: 

dX(/)  =  a(/,uo)d5(/)-bdZ(/).  (3) 

Likewise,  it  may  also  be  shown  that  the  CSD  matrix 
of  the  array  output  is  given  as 

dPM  =  E[dX{f)dX{f)»] 

=  a(/,uo)a(/,uo)"d.P,(/)-l-dPz(/)  (4) 

Then  the  power  spectrum  P,(/,  u)  of  the  classi¬ 
cal  delay-sum  DF  is  given  by  [IJ;  P*(/,u)  = 
a(/,u)^ WPx(/)Wa(/,  u),  where  W  is  the  diagonal 
matrix  of  the  weights.  As  a  special  case,  suppose 
the  sensor  noises  are  independent  and  identically  dis¬ 
tributed  so  that  Pzif)  is  a  scalar  multiple  of  the  iden¬ 
tity  matrix.  In  this  case,  Px(/,  uq)  is  the  maximum 
value  and  illustrates  the  idea  behind  DF:  that  the 
source  direction  can  be  determined  by  examining  the 
maximum  Px{f,u)  as  a  function  of  the  "look  vector” 
u  . 

As  pointed  out  by  Johnson  [10],  obtaining  the  CSD 
matrix  Px{f)  is  basic  to  this  classical  DF,  as  well  as 
the  MVDR  and  MUSIC  DF.  Therefore,  consistent  es¬ 
timation  of  the  CSD  matrix  is  a  key  point  in  these  DF 
methods  and  is  usually  justified  on  exphcit  WSS  as¬ 
sumptions  and  implicit  ergodic  assumptions.  Before 
we  consider  those  cases  where  the  array  output  {x(0} 
is  non-stationary,  the  next  section  is  devoted  to  showr 
ing  how  to  generalize  DF  via  the  power  spectrum  to 
DF  via  the  bispectrum. 

3  Bispectral  DF 

Higher-order  statistical  analysis  of  time  series  has 
been  undergoing  a  tremendous  development  in  the 
l2ist  decade.  Excellent  review  articles  of  this  area  are 
[13],  [14],  while  [19]  and  [18]  are  foundational  texts. 
In  addition,  Forster  and  Nikias  [3]  have  produced  a 
MVDR-like  DF  based  on  the  bispectrum  while  Porat 
and  Friedlander  [16]  DF  using  fourth-order  cumulants. 
In  this  section,  we  show  how  to  obtain  a  bispectral 
MUSIC  DF  and  a  mixed  bispectral  MUSIC  DF. 


Definition  3.1  [13]  If  {x(t)}  is  a  zero-mean,  third- 
order  stationary  random  process,  then  one  definition 
of  the  third-order  cumulant  is 


Cx(ri,r2)  =  jE[x(t  -t-ri)(gix(t  4-  r2)(8)x(t)], 

where  ®  is  the  Kronecker  product.  The  cross  bispec¬ 
tral  density  tensor  (CBT)  is: 

Px(/l./2)  = 

/OO  fOO 
•OO  J  ~oo 

The  generalization  of  equation  (3)  for  the  CSD  matrix 
of  the  array  output  {x(t)}  to  the  CBT  is  straight¬ 
forward  [19]: 

dSx(/i./2)  =  E[dX{h)®dX{h)®dX{h  +  /2)].  (5) 

For  DF  using  the  bispectrum,  the  classical  assump¬ 
tions  2.1  need  to  be  modified: 

Assumption  3.1  Single  source,  independent  Gaus¬ 
sian  noise: 

1.  the  plane-wave  signal  a(t)  «5  a  zero-mean,  third-order 
stationary  random  process. 

2.  the  noise  z(t)  is  a  zero-mean,  Gaussian,  1FS5  ran¬ 
dom  process. 

3.  the  signal  5(1)  and  the  sensor  noise  are  independent. 

The  bispectrum  of  a  Gaussian  WSS  is  zero  [19,  page 
36].  From  equation  (3),  equation  (5),  and  assump¬ 
tion  3.1,  it  can  be  shown  that 

dSx(/i./2)  = 

«(/i,uo)  ®  a(/2,Uo)  ®  a(/i  -|-  /2,uo)d5x(/i. /2). 

In  this  form,  the  third-order  signal  subspace  is  ob¬ 
vious.  For  computational  purposes,  we  can  use  the 
equivalent  form  for  three  tensors:  x®y®?~  (x® 
y)z^,  which  permits  the  application  of  the  SVD  to  de¬ 
termine  the  projections  for  the  third-order  signal  and 
noise  subspaces.  Moreover,  this  equivalence  combined 
with  equation  (5)  indicates  that  the  CBT  may  be  es¬ 
timated  by  averaging  Kronecker  products  of  the  DFT 
obtained  by  segmenting  the  data  into  time  blocks  [14]. 

Bispectral  MUSIC  performs  DF  by  projecting  a 
steering  vector  onto  the  third-order  noise  subspace  and 
calculating  the  reciprocal  of  the  norm  of  the  projec¬ 
tion.  Mixed  bispectral  MUSIC  performs  DF  by  pro¬ 
jecting  a  steering  vector  onto  the  second-order  noise 
subspace  fstandard  MUSIC)  followed  by  the  projec¬ 
tion  into  tne  third-order  noise  subspace  and  then  cal¬ 
culating  the  reciprocal  of  the  norm. 

The  following  simulation  provides  an  assessment  of 
these  DF  methods.  A  10-sensor  linear  array  is  laid 
along  the  z-axis.  Each  sensor  is  spaced  one-half  wave 
length  at  11  Hertz.  The  standard  right-hand  coordi¬ 
nate  system  is  assumed  so  that  broilside  to  the  ar¬ 
ray  measure  9  =  90®.  A  Gaussian  and  a  bilinear  sig¬ 
nal  were  present.  The  Gaussian  signal  was  generated 


by  AR  filtering  i.i.d.  A'lO,  1)  noise  (AR  coefficients; 
00=1,  ai=.9,  oj=.81).  This  signal  was  normalized  by 
setting  the  sample  variance  to  1.  It  arrived  on  the 
array  from  80®  azimuth  using  an  FFT  to  perform  the 
lag.  The  bilinear  signal  was  generated  using  the  bi¬ 
linear  model  of  Gabr  [5]:  St  -  1.3st_3Ct_3  -bf*,  where 
{ft}  is  i.i.d.  //’(0, 1)  and  then  normalized  by  setting 
the  sample  variance  to  1 .  It  arrived  on  the  array  from 
100®  azimuth  using  zin  FFT  to  perform  the  lag.  Fi¬ 
nally,  the  sensor  noise  was  modeled  as  i.i.d.  A/^O,  1). 
The  sample  rate  was  34.38  Hertz  and  a  total  of  16,384 
array  snapshots  were  collected. 

The  classical,  MUSIC,  bispectral  MUSIC  and 
mixed  bispectral  MUSIC  DF  were  applied  to  this  data. 
All  CSD  and  CBT  estimates  were  made  with  32- 
sample  block  averaging  assuming  /i=/2=11.5  Hertz. 
All  MUSIC-type  DF  assumed  two  signals  present.  The 
results  are  plotted  in  figure  1.  The  classical,  MUSIC, 
and  mixed  bispectral  MUSIC  indicate  both  sources 
with  the  primary  lobe  at  the  Gaussian  source  while  the 
bispectral  MUSIC  appears  to  discriminate  between 
the  Gaussian  and  the  bilinear  source. 


Figure  1;  Bispectral  DF  on  Two  Sources  —  Gaus¬ 
sian  source  at  80®  azimuth;  bilinear  source  at  100® 
azimuth;  Gaussian  sensor  noise.  The  solid  line  is  the 
classical  DF,  the  dashed  line  is  MUSIC,  the  dotted  line 
is  bispectral  MUSIC,  and  the  dash-dot  line  is  mixed 
bispectral  MUSIC. 

4  DF  via  the  Gabor  Transform 

This  section  shows  how  the  Gabor  transform  can 
be  used  as  a  pre-filter  for  removing  transient  signals 
corrupting  the  sensor  data.  Since  a  transient  can  be 
isolated  from  the  sensor  data,  it  follows  that  the  same 
processing  can  apply  to  transient  DF  [22,  Chapter 
4]  and  admit  more  advanced  transient  DF  based  on 
minimum-variance  or  subspace  methods.  For  brevity, 
we  focus  on  using  the  Gabor  transform  as  a  pre-filter. 

The  application  of  the  Gabor  transform  as  a  means 
of  transient  detection  was  beautifully  developed  by 
Friedlander  and  Porat  [4]  and  we  borrow  heavily  from 
that  paper.  Let  y(t)  be  a  given  continuous-time,  real¬ 
valued  signal.  Let  g{t)  be  a  fixed,  non-negative  func¬ 
tion  of  unit  area,  callea  the  window  function.  Then  the 


Gabor  representation  of  j/(<)  using  the  window  func¬ 
tion  g{t)  is 

fn=-^oon=:  — oo 

where  a  and  0  are  positive,  with  aP  <  1.  This  repre¬ 
sentation  shows  that  the  signal  is  to  be  expanded  into 
time-shifted  and  frequency-shifted  versions  of  the  win¬ 
dow.  Therefore,  the  application  will  dictate  the  choice 
of  the  window.  Indeed,  this  is  the  basic  point  [4]  for 
analyzing  transient  signals.  Since  they  consider  a  tran¬ 
sient  signal  as  a  signal  whose  duration  is  short  com¬ 
pared  to  the  observation  window,  their  basic  idea  is  to 
replace  the  original  Gaussian  window  of  Gabor  by  the 
one-sided  exponential  function:  g{t)  =  V^e~^*u(i) 
to  better  follow  transient  signals  with  a  sudden  “on” 
time  and  subsequent  exponential  decay.  Here  u(t)  is 
the  unit  step  function  and  “the  parameter  A  is  used 
to  control  the  effective  width  of  the  window”  [4]. 

The  use  of  the  Gabor  transform  with  this  window 
is  shown  in  the  following  DF  simulation.  A  10-sensor 
linear  array  is  laid  along  the  x-axis  and  each  sensor  is 
spaced  one-half  wave  length  at  11  Hertz  with  the  first 
sensor  located  at  the  origin.  An  11-Hertz  plane- wave 
signal  is  impinging  broadside  to  the  array  (zaimuth 
0  =  90®)  with  an  amplitude  of  1/10.  The  signal  is  cor¬ 
rupted  by  i.i.d.  Gaussian  noise  A  (0, 1)  at  each  sensor. 
In  addition,  five  one-sided  decaying  transients  were 
added  to  the  environment  (o=l,  /?=1,  A=l)  and  rele¬ 
vant  parzimeters  zire  displayed  in  Table  2:  A  total  of 


Table  2:  Parameters  for  Gabor  Transients 


index 

time  delay 
(sec) 

azimuth 

(deg) 

Irequcncy 

(Hertz) 

1 

13 

62 

12 

2 

23 

63 

12 

3 

21 

110 

10 

4 

5 

64 

12 

5 

10 

80 

13 

1,024  array  snapshots  were  collected  during  the  time 
interval  of  0  <  <  <  32  seconds  at  a  sampling  rate  of  32 
Hertz.  Figure  2  plots  the  11-Hertz  signal,  the  Gaus¬ 
sian  noise,  euid  the  five  transients  at  the  first  sensor. 
Figure  3  is  a  plot  of  the  discrete  Gabor  transform  of 
the  output  of  first  sensor.  The  time  axis  starts  at  time 
t  =  0  in  the  leftmost  corner  and  moves  in  one-second 
increments  to  the  uppermost  corner.  The  frequency 
axis  starts  at  /  =  0  in  the  leftmost  corner  and  moves 
in  one-Hertz  increments  to  the  lowest  corner.  The  five 
transients  are  clearly  displayed. 

The  effect  of  the  transients  on  the  MVDR  DF  is 
shown  in  Figure  4.  MVDR  used  11  Hertz  in  an  at¬ 
tempt  to  find  the  weak  11-Hertz  signal  but  the  tran¬ 
sients  swamp  the  output.  The  transients  were  removed 
by  tzdcing  the  Gabor  transform  of  each  sensor’s  output, 
finding  those  coefficients  exceeding  three  standard  de¬ 
viations  of  the  Gabor  transform,  and  then  subtracting 


0| 


the  corresponding  Gabor  transients  (weighted  by  its 
Gabor  coefEcientJ.  Thus,  the  Gabor  transform  func¬ 
tioned  as  a  pre-niter  and  the  resulting  array  output 
was  passed  to  the  DF.  With  the  transients  removed, 
the  MVDR  DF  located  the  11-Hertz  signal. 

»r - - - ■ - ^ ^ - - - ■ - - - 1 


6 


- - - - - - - ^ - - - . - 1 

0  5  10  15  20  25  30  35 

bac  (s«c) 


Figure  2:  Time  Series  of  the  First  Sensor  —  The 
11-Hertz  signal  at  1/10  amplitude,  additive  ,^(0,1) 
noise,  and  the  five  Gabor  transients. 


Figure  3;  Gabor  Transform  of  the  Time  Series 
of  the  First  Sensor  —  The  11-Hertz  signal  at  1/10 
amplitude,  additive  ^/{0, 1)  noise,  and  the  five  Gabor 
transients. 


5  DF  via  the  WVD 

This  section  shows  how  one  preliminary  estimate 
of  the  Wigner-Ville  Distribution  (WVD)  of  the  ar¬ 
ray  output  can  be  used  to  DF  a  source  in  both  time 
and  frequency.  This  approach  is  different  from  that  of 
Swindlehurst  and  Kailath  [21]  in  which  a  spatial  WVD 
was  used  to  handle  the  near-field  effects  and  produced 
both  bearing  and  range  estimates.  Instead,  we  show 
how  to  use  the  WVD  to  DF  on  a  source  in  the  far  field, 
provided  the  signal  admits  a  WVD  sufficiently  differ¬ 
ent  from  the  background  noise  at  a  particular  time 
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Figure  4:  MVDR  using  a  Gabor  PrefUter  —  An 
11-Hertz  simal  at  1/10  amplitude  from  90®  azimuth, 
additive  A^(0,1)  noise,  and  the  five  Gabor  transients. 
The  CSD  was  estimated  at  1 1  Hertz  by  averaging  32- 
sample  blocks  of  the  1,024  array  snapshots.  The  solid 
line  is  MVDR  without  the  Gabor  prefilter  and  the 
dashed  is  MVDR  with  the  Gabor  prefilter. 


and  frequency.  As  such,  this  approach  could  be  used 
to  DF  on  transients. 

Given  two  analytic  signals  x(t)  and  y(t),  their  cross 
Wigner-Ville  Distribution  (XWVD)  is  [2]: 


VF(x,y,t, 


c-'2'^i(t-l-T/2)y(t~r/2)dT. 


Consider  the  case  of  the  array  output  when  only  one 
analytic  signal  is  present  from  the  direction  uq.  Then 
x(t)  =  [s(t  -f  Tm(uo))].  The  XWVD  of  the  mth  and 
m'th  sensor  outputs  is: 

W(r„,  t,  /)  = 

X  iy(5,5;<-|-r,„(uo)/2-f-r„.(uo)/2,/). 

If  the  XWVD  is  smoothed  in  time,  then  three  obser¬ 
vations  can  be  made:  First,  if  the  time  lags  are  rel¬ 
atively  small  with  respect  to  the  smoothing  window, 
then  the  phase  term  is  not  affected  while  the  time- 
shifted  WVD  of  the  signal  may  be  approximated  by 
the  WVD  of  the  signal.  In  matrix  form,  we  get  the 
estimate  of  the  WVD  of  the  array: 


W(x,x,i,f)  w  a(/, uo)a(/,Uo)"W(s,s;t,/), 


where  W  denotes  the  smoothed  WVD  of  the  array 
output.  Second,  if  noise  is  present  and  distributes 
in  a  uniform  and  zero-mean  fashion  over  the  time- 
frequency  plane,  then  the  smoothing  should  have  the 
effect  of  providing  a  biased  estimate  of  the  WVD  of 
the  array  but  with  a  lower  variance  than  the  “raw” 
WVD.  Third,  the  form  of  the  estimate  permits  the 
application  of  subspace  methods. 

The  following  simulation  supports  these  claims. 
The  10-sensor  linear  array  spaced  at  1/2  wavelength 


(/  =  11  Hertz)  is  receiving  chirp  signals  from  two 
sources:  one  source  is  broadside  to  the  array  (azimuth 
6  =  PO®).  and  transmits  a  four-second  ciirp.  The 
chirp  starts  at  time  t  =  0  at  16  Hertz  and  sweeps 
down  at  4  Hz/sec.  The  second  source  is  located  at 
6  =  60®  azimuth  and  transmits  a  chirp  sweeping  up  at 
4  Hz/sec  starting  from  0  Hertz.  Sensor  noise  was  mod¬ 
eled  as  i.i.d  .A/’(0,  1).  A  toted  of  128  array  snapshots 
were  collected  at  a  32  Hertz  sample  rate.  Figure  5 
plots  conventional  DF  at  11  Hertz,  the  conventional 
WVD  DF,  and  the  MUSIC  WVD  (one  signal).  Both 
WVD  DF  used  time  t  =  1.25  seconds  and  /  =  11 
Hertz.  The  conventional  DF  estimated  the  CSD  ma¬ 
trix  using  32-sample  time  blocks.  The  XWVD  was 
estimated  an  FFT  size  of  64  and  rectangular  8-point 
time-smewthing.  The  figure  shows  that  both  the  WVD 
DF  distinguish  the  two  sources  while  the  conventional 
DF  also  shows  that  there  are  two  signals  at  11  Hertz 
at  60®  and  90®  azimuth. 


Figure  5;  WVD  DF  on  Two  Chirps  —  One  chirp 
at  90®  azimuth  sweeping  down  from  16  Hertz  at  -4 
Hz/sec.  The  second  chirp  at  60®  azimuth  sweeping 
up  from  0  Hertz  at  -)-4  Hz/sec.  The  solid  line  is  the 
classical  DF,  the  dashed  line  is  classical  WVD  DF,  and 
the  dotted  line  is  the  MUSIC  WVD  DF  (one  signal). 


6  DF  on  Cyclostationary  Signals 

The  fundamental  stumbling  block  in  treating  non¬ 
stationary  time  series  is  the  lack  of  either  a  sample 
ensemble  (which  would  permit  a  consistent  estimation 
of  the  CSD  matrix),  or  an  ergodic  theorem  (which  can 
produce  a  consistent  estimate  from  a  sufficiently  long 
realization  of  the  time  series).  While  there  is  a  large 
literature  of  attacks  on  non-stationary  time  series  (see 
Priestley  fl7]  for  an  accessible  survey),  this  section 
will  consider  harmonizable,  cyclostationary  time  se¬ 
ries.  These  assumptions  will  permit  the  estimation  of 
the  CSD  matrix  and  readily  apply  to  DF. 

Cyclostationary  time  series  have  been  undergoing 
an  extensive  development  in  the  last  decade,  (the 
April  1991  issue  of  the  IEEE  Signal  Processing  Mag- 
azine  is  devoted  to  cyclostationary  processing).  In 
particular,  we  single  out  the  excellent  Ph.D  thesis 


of  Schell  [20]  which  exploits  cyclostationarity  to  per¬ 
form  DF  using  subspace  methods  and  has  performed  a 
Cramer-Rao  lower  bound  analysis.  Another  approach 
to  DF  can  be  made  using  spectral  estimators  obteiined 
by  Hurd  [7],  [8]  and  is  the  subject  of  this  section. 

Definition  6.1  [15,  page  226]  A  random  process 
(z^t)}  is  called  cyclostationary  of  period  T  >  0  for 
which  the  mean  and  covariance  functions  are  periodic: 
for  all  times  t,  ti,  and  tj,  there  holds  E\x{t  -f  T)]  = 
F^[i(t)]  and  Rx(ti ,  t2)  =  Rx{ti  -I-  T,  <2  -t-  T). 

Definition  C.2  [12]  A  random  process  {a:(t)}  is  called 
(strongly)  harmonizable  provided  it  can  be  repre¬ 
sented  in  quadratic  mean  for  every  time  t  by  the  inte¬ 
gral 

/oo 

e‘^/‘dX(/) 

■OO 

where  X  is  a  random  measure  for  which  the  spectral 

measure  defined  by  Px{A  x  B)  =  E[X{A)X{B)]  is  of 
bounded  variation. 

Hurd  [7]  points  out  that  when  {*(0)  is  both  cy¬ 
clostationary  and  harmonizable  then  the  spectrum 
Pxifi^fi)  is  confined  to  diagonal  lines  in  the  fi  x  /2 
plane  which  are  p2irallel  to  the  main  diagonal  and 
spaced  l/T  Hertz  apart.  In  particular,  stationary 
noise  will  be  confined  to  the  main  diagonal  (/i=/2) 
and  DF  using  fx  ^  /2  should  permit  noise  suppres¬ 
sion.  The  DF  assumptions  are  as  follows: 

Assumption  6.1  Single  Cyclostationary  source, 
uncorrelated  noise: 

1.  the  plane  wave  signal  5(t)  tr  a  strongly  harmonizable, 
cyclostationary  random  process. 

Z.  the  noise  z(t)  is  zero-mean,  WSS. 

3.  the  signal  s{t)  and  the  sensor  noise  are  uncorrelated. 

Since  the  noise  vector  z(t)  is  WSS  and  signal  s(t)  is 
strongly  harmonizable,  then  the  array  Then  the  array 
output  x(f)  is  harmonizable  with  a  Fourier  measure 
given  by  dX{f)  =  a(/,  uo)d5f/)  -f  dZ(/),  where  only 
Z(/)  is  orthogonally  scatterea.  From  this  representa¬ 
tion,  the  CSD  matrix  of  the  array  hee  the  form: 

dPx{fi,f2)  =  E[dX{fx)dX{f2)"] 

=  <fP.(/i,/2)a(/i,uo)a(/2,uo)"  -bdPn(/l,/2). 

Note  the  following:  the  noise  CSD  matrix  is  supported 
only  on  the  line  of  zero  frequency  difference;  the  exten¬ 
sion  to  subspace  algorithms  appears  straight-forward; 
the  CSD  matrix  can  be  estimated  by  averaging  outer 
products  of  FFT  vectors  of  the  time  segments  of  the 
data  if  the  period  T  is  known. 

The  following  simulation  examines  these  claims. 
A  10-sensor  linear  array  spaced  at  1/2  wavelength 
(/  =  11  Hertz)  is  receiving  two  signals:  the  first 
is  a  cyclostationary  signal  generated  by  amplitude- 
modulating  Af(0, 1)  noise  at  4  Hertz.  This  AM  signal 
is  normalized  by  setting  the  sample  standard  devia¬ 
tion  to  one  and  arrives  broadside  to  the  array  (90® 


azimuth).  The  second  signal  is  an  11-Hertz  complex 
exponential  arriving  at  75°  azimuth.  Sensor  noise  is 
i.i.d.  Af(0, 1).  A  total  of  256  array  snapshots  were 
collected  at  a  128  Hertz  sample  rate.  The  classical 
DF  estimated  the  CSD  matrix  by  averaging  32-point 
blocks  at  11  Hertz.  The  cyclostationary  DF  worked  by 
picking  a  look  direction  and  then  feeding  the  output 
of  the  conventional  delay-sum  beamformer  into  a  esti¬ 
mator  of  the  two-dimensional  spectrum  ixC/ii/z)  at 
the  point  (/i./z)  =  (H.  15)  Hertz.  The  spectral  esti¬ 
mator  performed  diagonal  smoothing  by  averaging  11 
diagonal  bins  spaced  0.5  Hertz  apart  centered  about 
at  (11,15).  The  results  are  displayed  in  figure  6.  The 
classical  DF  localizes  the  11-Hertz  signal  while  the  cy¬ 
clostationary  DF  locates  the  AM  signal. 


Figure  6;  DF  on  Cyclostationary  Signals  — 
Gaussian  AM  signal  at  90°  azimuth;  11  Hertz  signal 
at  75°  azimuth.  The  solid  line  is  tiie  classical  DF;  the 
dashed  line  is  the  cyclostationary  DF. 

7  DF  on  a  Chaotic  Signal 

In  the  past  decade,  there  has  been  a  surge  of  pub¬ 
lications  in  ch2U3tic  systems.  Of  particular  interest 
for  sonar  applications,  we  single  out  the  work  of  K. 
Karagiannis  [11]  which  observed  chaotic  behavior  in 
the  noise  generated  by  rattling  gearboxes.  One  basic 
model  of  cheiotic  behavior  is  the  conventional  logistic 
map  (CLM)  [23]: 

x(t-(-  1)  =  ;jx(t)(l -z(t)). 

One  stochastic  generalization  of  the  CLM  is  the  Model 
2  of  [6]:  for  t  =  1,  2,  ... 

x(( -Kl)  -  k(a,6)x(0“(l  -  x(0)‘«(«  +  1)', 

where  a,  b,  and  c  are  non-negative,  {«(<)}  are  inde¬ 
pendent  and  identically  distributed  uniformly  on  the 
interval  [0, 1],  k{a,b)  =  (a -i- 6)“"''‘/(a“6*),  and  the  ini¬ 
tial  value  is  set  as  x(0)  =  u(0)°.  Estimators  of  a,  b,  c 
were  obtained  and  it  was  observed  that  Model  2  pro¬ 
vided  a  good  fit  to  time  series  generated  by  the  CLM. 

This  section  demonstrates  that  the  Model  2  esti¬ 
mators  of  a,  b,  and  c  can  apply  to  source  localization. 


The  basic  idea  is  to  assume  that  a  signal  s(t)  is  of  the 
form  of  the  Model  2.  A  look  direction  is  picked  and  the 
time  series  resulting  from  the  conventional  delay-sum 
beamforming  is  fed  into  the  estimation  scheme.  The 
one-step  mean-squared  prediction  error  is  calculated 
and  its  reciprocal  is  plotted  as  a  function  of  the  look 
direction.  If  the  look  direction  matches  the  signal  di¬ 
rection,  then  a  good  fit  of  the  signal  is  obtained  which 
results  in  a  small  prediction  error.  Otherwise,  a  poor 
fit  occurs  which  results  in  a  large  prediction  error.  We 
compare  this  FractaJ  DF  with  the  conventional  DF 
and  the  MUSIC  DF  in  the  following  simulations. 

We  use  a  10-sensor  linear  array  with  1/2  wavelength 
spacing  (/  =  11  Hertz).  One  signal  is  placed  broad¬ 
side  to  the  array  (azimuth  0  =  90°)  and  is  generated 
from  the  CLM  (/i  =  4,  s(0)  =  1/3).  An  second  signal 
was  generated  from  the  CLM  and  arrives  on  the  array 
from  95®  aizimuth.  The  arrival  was  simulated  using  the 
Fourier  method:  the  FFT  of  the  time  series  was  mul¬ 
tiplied  by  the  appropriate  phase  for  each  sensor  and 
the  inverse  FFT  determined  the  signal  at  each  sensor. 
Sensor  noise  was  modeled  as  i.i.d  zero  mean  Gaussian 
with  the  standard  deviation  set  so  the  SNR=1.  A 
total  of  17=128  array  snapshots  were  collected  at  a 
100  Hertz  sample  rate.  The  Classical  and  MUSIC  DF 
were  operated  at  11  Hertz  and  the  CSD  matrix  was 
estimated  using  32-point  block-averaging.  MUSIC  as¬ 
sumed  two  sources.  All  DFs  locate  the  source  at  90° 
azimuth  though  the  beamwidth  of  the  fractal  DF  is 
much  narrower  than  the  other  DFs.  The  results  are 
shown  in  figure  7.  The  Conventional  and  MUSIC  DFs 
could  not  resolve  this  sources,  but  the  Fractal  DF  is 
correctly  indicating  the  signal  directions  at  90°  and 
95°  azimuth. 


Figure  7:  DF  on  a  Two  Chaotic  Signals  —  Two 
CLM  signals  at  90°  and  95°  azimuth.  The  solid  line 
is  the  classical  DF,  the  dashed  line  is  MUSIC,  and  the 
dotted  line  is  the  Fractal  DF. 


8  Conclusions 

The  preceding  sections  have  shown  that  Johnson’s 
[10]  approach  to  bearing  estimation  as  form  of  spectral 
estunation  readily  generalizes  whenever  a  statistically 
ronsisf<*nt  estimator  of  a  signal  attribute  c^u  be  found. 
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